Cargamos los arhivos contenidos en “infoClin_all_dataset.RData” y “CRC_stageII.RData” con los contenidos clínicos y de CRC.
# Cargamos el archivo infoClin_all_dataset.RData 16/03/2024
dataClin <- load("data/infoClin_all_dataset.RData")
# Cargamos el archivo CRC_stageII.RData 19/03/2024 dataCRC <-
# load('data/CRC_stageII.RData')
# Cargamos el archivo exp_mat.RData 02/04/2024
dataCRC <- load("data/exp_mat.RData")
GSE Files
Platform GPL570: [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array
Fusionamos los datos en un único dataframe al que le añadimos la columna grupo con los datos stage y msi.
Factorizamos las variables categóricas y ordenamos por ID tanto los datos clínicos como los de la matriz de expresión.
Creamos un DataFrame para asignar cada ID a su estudio clínico. A cada estudio le asignamos un color para poder identificarlo fácilmente en los estudios gráficos.
# 11/04/2023 Función para preparar el data frame con el estudio por ID. También
# asignaremos un color a cada estudio. Utilizamos el paquete 'dplyr'
prepare_df <- function(df, study_name, color) {
df %>%
select(ID) %>%
mutate(study = study_name, color = color)
}
# Lista de estudios con sus colores asociados
studies_info <- list(clx = "red", GSE39582 = "blue", GSE14333 = "green", GSE13294 = "orange",
GSE17536 = "yellow", TCGA = "violet")
# Creamos una lista de de los estudios utilizando la función prepare_df
df_studies_list <- lapply(names(studies_info), function(x) {
prepare_df(get(paste0(x, "_infoClin")), x, studies_info[[x]])
})
# Juntamos los df's de la lista se studies en uno solo y lo ordenamos por 'ID'
df_studies <- bind_rows(df_studies_list) %>%
arrange(ID)
# Fusionamos y ordenamos todos los datos clínicos generando 'df_dataClin'
df_dataClin <- bind_rows(lapply(names(studies_info), function(x) get(paste0(x, "_infoClin")))) %>%
arrange(ID)
# Cramos una nueva columna 'grupo' a los 'df_dataClin' con los datos de msi y
# stage
df_dataClin$grupo <- with(df_dataClin, paste(stage, msi_imputed, sep = "_"))
# Cramos una nueva columna 'cms_stage' a los 'df_dataClin' con los datos de cms
# y stage 19/04/2024
df_dataClin$cms_stage <- with(df_dataClin, ifelse(is.na(cms), NA, paste(stage, cms,
sep = "_")))
# Convertimos a factores
df_dataClin <- df_dataClin %>%
mutate(across(c(msi_imputed, cms, stage, grupo, cms_stage), factor))
# Añadimos la columna 'study' a df_dataClin. Al final no añadimos 'color'
df_dataClin <- left_join(df_dataClin, select(df_studies, ID, study), by = "ID")
Pasamos la matriz de expresión a data_frame y la ordenamos por “ID”
# Pasamos la matriz de expresión a data frame 26/03/2024
df_exmat <- as.data.frame(ex)
# Ordenamos la matriz de expresión por nombre de columnas 26/03/2024
df_exmat <- df_exmat[, order(names(df_exmat))]
# Transponemos la matriz de expresión: 22/04/2025
df_exmat_t <- as.data.frame(t(df_exmat))
# Creamos un df con los datos clínicos y los de expresión génica (transpuesta)
df_dataClin_exmat_t <- cbind(df_dataClin, df_exmat_t)
Resumen exploratorio de los datos:
# Hacemos un Summary con el paquete skimr 18/03/2024
skim(df_dataClin)
| Name | df_dataClin |
| Number of rows | 1079 |
| Number of columns | 7 |
| _______________________ | |
| Column type frequency: | |
| character | 2 |
| factor | 5 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| ID | 0 | 1 | 7 | 12 | 0 | 1079 | 0 |
| study | 0 | 1 | 3 | 8 | 0 | 6 | 0 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| msi_imputed | 0 | 1.00 | FALSE | 2 | MSS: 840, MSI: 239 |
| cms | 158 | 0.85 | FALSE | 4 | CMS: 329, CMS: 262, CMS: 194, CMS: 136 |
| stage | 0 | 1.00 | FALSE | 2 | II: 684, III: 395 |
| grupo | 0 | 1.00 | FALSE | 4 | II_: 511, III: 329, II_: 173, III: 66 |
| cms_stage | 158 | 0.85 | FALSE | 8 | II_: 210, II_: 151, II_: 124, III: 119 |
Estadísticas y cálculos para la generación de los gráficos de datos clínicos:
# Estadísticas desglosadas 18/03/2024 Frecuencias
msi_frecuencias <- table(df_dataClin$msi_imputed)
cms_frecuencias <- table(df_dataClin$cms)
stage_frecuencias <- table(df_dataClin$stage)
# Valores NA
na_msi <- sum(is.na(df_dataClin$msi_imputed))
na_cms <- sum(is.na(df_dataClin$cms))
na_stage <- sum(is.na(df_dataClin$stage))
# Contamos cuantos MSI y MSS tenemos por cada estadio, II y III 19/03/2024
sII_MSI <- sum(df_dataClin$msi_imputed == "MSI" & df_dataClin$stage == "II")
sII_MSS <- sum(df_dataClin$msi_imputed == "MSS" & df_dataClin$stage == "II")
sIII_MSI <- sum(df_dataClin$msi_imputed == "MSI" & df_dataClin$stage == "III")
sIII_MSS <- sum(df_dataClin$msi_imputed == "MSS" & df_dataClin$stage == "III")
# Crear un dataframe con estos valores
stage_MSI_MSS_df <- data.frame(Stage_MSI_MSS = c("II_MSI", "II_MSS", "III_MSI", "III_MSS"),
Freq = c(sII_MSI, sII_MSS, sIII_MSI, sIII_MSS))
# Contamos cuantos MSI y MSS tenemos por cada estadio, II y III 19/04/2024
sII_CMS1 <- sum(df_dataClin$cms_stage == "II_CMS1", na.rm = TRUE)
sII_CMS2 <- sum(df_dataClin$cms_stage == "II_CMS2", na.rm = TRUE)
sII_CMS3 <- sum(df_dataClin$cms_stage == "II_CMS3", na.rm = TRUE)
sII_CMS4 <- sum(df_dataClin$cms_stage == "II_CMS4", na.rm = TRUE)
sIII_CMS1 <- sum(df_dataClin$cms_stage == "III_CMS1", na.rm = TRUE)
sIII_CMS2 <- sum(df_dataClin$cms_stage == "III_CMS2", na.rm = TRUE)
sIII_CMS3 <- sum(df_dataClin$cms_stage == "III_CMS3", na.rm = TRUE)
sIII_CMS4 <- sum(df_dataClin$cms_stage == "III_CMS4", na.rm = TRUE)
# Crear un dataframe con estos valores
stage_CMS_df <- data.frame(Stage_CMS = c("II_CMS1", "II_CMS2", "II_CMS3", "II_CMS4",
"III_CMS1", "III_CMS2", "III_CMS3", "III_CMS4"), Freq = c(sII_CMS1, sII_CMS2,
sII_CMS3, sII_CMS4, sIII_CMS1, sIII_CMS2, sIII_CMS3, sIII_CMS4))
# Pasamos a DataFrame stage
df_stage_frecuencias <- as.data.frame(stage_frecuencias)
df_stage_frecuencias
Var1 Freq
1 II 684
2 III 395
# Añadimos los porcentajes
# Calculamos el total de frecuencias
total_frecuencias <- sum(df_stage_frecuencias$Freq)
# Añadimos una columna de porcentajes al dataframe
df_stage_frecuencias$Porcentaje <- (df_stage_frecuencias$Freq / total_frecuencias) * 100
# Diagrama de Barras de ESTADIO
p1 <- ggplot(df_stage_frecuencias, aes(x = Var1, y = Freq)) +
geom_bar(stat = "identity", fill = c("salmon", "seashell3")) +
geom_text(aes(label = Freq), vjust = -0.5, size = 3.5) + # Frecuencias arriba de cada barra
# Porcentajes en la mitad + o -
geom_text(aes(label = paste0(round(Porcentaje, 1), "%")), vjust = 5, size = 3.5) +
labs(title = "STAGE", x = "", y = "Frecuencia") +
theme_minimal()
p1
# Pasamos a DataFrame msi
df_msi_frecuencias <- as.data.frame(msi_frecuencias)
df_msi_frecuencias
Var1 Freq
1 MSI 239
2 MSS 840
# Añadimos los porcentajes
# Calculamos el total de frecuencias
total_frecuencias <- sum(df_msi_frecuencias$Freq)
# Añadimos una columna de porcentajes al dataframe
df_msi_frecuencias$Porcentaje <- (df_msi_frecuencias$Freq / total_frecuencias) * 100
# Diagrama de Barras de MSI/MSS
p2 <- ggplot(df_msi_frecuencias, aes(x = Var1, y = Freq)) +
geom_bar(stat = "identity", fill = c("salmon", "seashell3")) +
geom_text(aes(label = Freq), vjust = -0.5, size = 3.5) + # Frecuencias arriba de cada barra
# Porcentajes en la mitad + o -
geom_text(aes(label = paste0(round(Porcentaje, 1), "%")), vjust = 5, size = 3.5) +
labs(title = "MSI/MSS", x = "", y = "Frecuencia") +
theme_minimal()
p2
# Stage vs MSI/MSS DataFrame
stage_MSI_MSS_df
Stage_MSI_MSS Freq
1 II_MSI 173
2 II_MSS 511
3 III_MSI 66
4 III_MSS 329
# Añadimos los porcentajes
# Calculamos el total de frecuencias
total_frecuencias <- sum(stage_MSI_MSS_df$Freq)
# Añadimos una columna de porcentajes al dataframe
stage_MSI_MSS_df$Porcentaje <- (stage_MSI_MSS_df$Freq / total_frecuencias) * 100
# Diagrama de Barras de Stage vs MSI/MSS
p3 <- ggplot(stage_MSI_MSS_df, aes(x = Stage_MSI_MSS, y = Freq)) +
geom_bar(stat = "identity", fill = c("salmon", "seashell3",
"turquoise3", "slategray")) +
geom_text(aes(label = Freq), vjust = -0.5, size = 3.5) + # Frecuencias arriba de cada barra
# Porcentajes en la mitad + o -
geom_text(aes(label = paste0(round(Porcentaje, 1), "%")), vjust = 2, size = 3.5) +
labs(title = "STAGE & MSI/MSS", x = "", y = "Frecuencia") +
theme_minimal()
p3
# Pasamos a DataFrame cms
df_cms_frecuencias <- as.data.frame(cms_frecuencias)
df_cms_frecuencias
Var1 Freq
1 CMS1 194
2 CMS2 329
3 CMS3 136
4 CMS4 262
# Añadimos los porcentajes
# Calculamos el total de frecuencias
total_frecuencias <- sum(df_cms_frecuencias$Freq)
# Añadimos una columna de porcentajes al dataframe
df_cms_frecuencias$Porcentaje <- (df_cms_frecuencias$Freq / total_frecuencias) * 100
# Diagrama de Barras de CMS
p4 <- ggplot(df_cms_frecuencias, aes(x = Var1, y = Freq, fill = Var1)) +
geom_bar(stat = "identity") +
geom_text(aes(label = Freq), vjust = -0.3, size = 3.5) + # Frecuencias arriba de cada barra
# Porcentajes en la mitad + o -
geom_text(aes(label = paste0(round(Porcentaje, 1), "%")), vjust = 5, size = 3.5) +
scale_fill_manual(values = c("salmon", "seashell3", "turquoise3", "slategray")) +
labs(title = "CMS", x = "", y = "Frecuencia") +
theme_minimal() +
theme(legend.position = "none") # Omitir la leyenda
p4
# Stage vs MSI/MSS DataFrame
stage_CMS_df
Stage_CMS Freq
1 II_CMS1 124
2 II_CMS2 210
3 II_CMS3 96
4 II_CMS4 151
5 III_CMS1 70
6 III_CMS2 119
7 III_CMS3 40
8 III_CMS4 111
# Añadimos los porcentajes
# Calculamos el total de frecuencias
total_frecuencias <- sum(stage_CMS_df$Freq)
# Añadimos una columna de porcentajes al dataframe
stage_CMS_df$Porcentaje <- (stage_CMS_df$Freq / total_frecuencias) * 100
# Diagrama de Barras de Stage vs MSI/MSS
p5 <- ggplot(stage_CMS_df, aes(x = Stage_CMS, y = Freq)) +
geom_bar(stat = "identity", fill = c("salmon", "seashell3",
"turquoise3", "slategray",
"salmon", "seashell3",
"turquoise3", "slategray")) +
geom_text(aes(label = Freq), vjust = -0.5, size = 3.5) + # Frecuencias arriba de cada barra
# Porcentajes en la mitad + o -
geom_text(aes(label = paste0(round(Porcentaje, 1), "%")), vjust = 2, size = 3.5) +
labs(title = "STAGE & CMS", x = "", y = "Frecuencia") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))# Inclinar las etiquetas del eje x
p5
# Múltiple Layout con patchwork En un grid 3x2
p2 + p3 + p4 + p5 + p1 + plot_layout(ncol = 2, nrow = 3)
Para el Análisis de Componentes Principales (PCA), normalizamos y
escalamos (estandarizamos) la matriz de expresión. En nuestro caso, los
datos ya están normalizamos (en principio con \(log_2\)) por lo que tan solo los
escalamos.
La estandarización se realiza restando la media de cada variable y
dividiendo el resultado por la desviación estándar de esa variable.
\[
Z = \frac{x-\mu}{\sigma}
\] Utilizamos la función scale para escalar. A
continuación realizamos el PCA con la transpuesta de la matriz de
expresión dado que queremos un análisis de las muestras en función de
los genes y no al revés.
IMP: al hacer el escalado antes de la transposición
estamos estandarizando por genes y no por muestras.
# Escalamos/normalizamos la matriz de expresión. También centra. 28/04/2024 Al
# hacer el escalado antes de la transposición estamos escalando por genes y no
# por muestras.
df_exmat_scaled <- scale(df_exmat)
# Hacemos el PCA. Hay que hacer la transpuesta de la matriz de expresión (las
# columnas representan genes y las filas representen muestras), así los
# componentes principales reflejan la variabilidad entre muestras. 28/04/2024
# Na haría falta center = TRUE ya que al escalar ya centralizamos.
pca_result <- prcomp(t(df_exmat_scaled), center = TRUE, scale. = FALSE)
RESUMEN PCA
# Obtenemos los valores propios (la varianza explicada por cada componente)
varianza <- pca_result$sdev^2
# Calculamos la proporción de varianza explicada por cada componente
prop_varianza <- varianza/sum(varianza)
# Calculamos la proporción acumulativa de varianza explicada
prop_acumulada <- cumsum(prop_varianza)
# Crear un dataframe con esta información para las primeras 10 PC
resumen_pca <- data.frame(SD = pca_result$sdev[1:10], Varianza = prop_varianza[1:10],
VarAcumulada = prop_acumulada[1:10])
print(resumen_pca)
SD Varianza VarAcumulada
1 13.721545 0.08197744 0.08197744
2 11.124090 0.05387873 0.13585617
3 10.811433 0.05089263 0.18674880
4 9.600035 0.04012675 0.22687555
5 6.928531 0.02090117 0.24777672
6 6.605919 0.01900005 0.26677677
7 6.062676 0.01600358 0.28278035
8 5.642600 0.01386267 0.29664302
9 5.370890 0.01255975 0.30920276
10 5.178365 0.01167545 0.32087821
La baja varianza explicada por los primeros dos componentes, PC1
(8,2%) y PC2 (5.4%), sugiere que los datos tienen una distribución de
varianza muy dispersa y que no hay unas pocas direcciones principales de
variabilidad que dominen.
Esto es coherente con el hecho de que estamos tratando con datos
biológicos, en concreto de CRC, que suelen ser de alta dimensionalidad
con múltiples factores contribuyendo a la variabilidad general. Es
decir, no hay patrones dominantes de variabilidad que puedan ser
fácilmente capturados por los primeros componentes principales.
CARGAS DE GENES EN PCA
Para entender mejor qué genes están contribuyendo a la separación observada, podemos analizar las cargas de los componentes principales para ver qué genes tienen mayores pesos en los componentes que diferencian los grupos.
# Extracción de cargas de componentes principales
cargas_pca <- pca_result$rotation
# Mostramos las cargas de los primeros dos componentes
primer_componente <- cargas_pca[, 1]
segundo_componente <- cargas_pca[, 2]
PC1
# Extraemos los top 10 genes y sus cargas de PC1
top_genes_pc1 <- head(sort(primer_componente, decreasing = TRUE), 10)
# Convertimos a data frame para visualización de PC1
top_genes_pc1_df <- data.frame(Gen = names(top_genes_pc1), Carga = top_genes_pc1)
# Ordenamos el dataframe por carga de manera descendente
top_genes_pc1_df <- top_genes_pc1_df[order(top_genes_pc1_df$Carga, decreasing = TRUE),
]
# Ajustamos los factores para que el orden de los genes sea según las cargas
top_genes_pc1_df$Gen <- factor(top_genes_pc1_df$Gen, levels = top_genes_pc1_df$Gen)
# Creamos el gráfico de barras
plot_pc1 <- ggplot(top_genes_pc1_df, aes(x = Gen, y = Carga, fill = Gen)) + geom_bar(stat = "identity",
position = "dodge") + theme_minimal() + theme(axis.text.x = element_text(angle = 90,
hjust = 1, vjust = 0.5)) + labs(title = "Top 10 Cargas de Genes en PC1", x = "Gen",
y = "Carga") + scale_color_gradient(low = "peachpuff", high = "peachpuff4")
# Mostramos los top 10 genes de PC1
print(top_genes_pc1)
SFRP2 GAS1 COL10A1 CCL18 CYP1B1 SPP1 GPNMB
0.06248246 0.04743630 0.04494118 0.04386073 0.04165178 0.03902618 0.03896303
GREM1 CCDC80 POSTN
0.03860923 0.03848662 0.03844061
PC2
# Extraemos los top 10 genes y sus cargas de PC2
top_genes_pc2 <- head(sort(segundo_componente, decreasing = TRUE), 10)
# Convertimos a data frame para visualización de PC2
top_genes_pc2_df <- data.frame(Gen = names(top_genes_pc2), Carga = top_genes_pc2)
# Ordenamos el dataframe por carga de manera descendente
top_genes_pc2_df <- top_genes_pc2_df[order(top_genes_pc2_df$Carga, decreasing = TRUE),
]
# Ajustamos los factores para que el orden de los genes sea según las cargas
top_genes_pc2_df$Gen <- factor(top_genes_pc2_df$Gen, levels = top_genes_pc2_df$Gen)
# Creamos el gráfico de barras
plot_pc2 <- ggplot(top_genes_pc2_df, aes(x = Gen, y = Carga, fill = Gen)) + geom_bar(stat = "identity",
position = "dodge") + theme_minimal() + theme(axis.text.x = element_text(angle = 90,
hjust = 1, vjust = 0.5)) + labs(title = "Top 10 Cargas de Genes en PC2", x = "Gen",
y = "Carga") + scale_color_gradient(low = "peachpuff", high = "peachpuff4")
# Mostramos los top 10 genes de PC1
print(top_genes_pc2)
REG4 ZIC2 REG1A TRIM7 CTSE TCN1 AGR3
0.07056229 0.05049736 0.04890627 0.04815804 0.04715623 0.04230700 0.04170686
RAB27B SDR16C5 PLA2G2A
0.04125286 0.04028539 0.03911398
# Mostramos gráficamente las cargas del top 10 tanto de PC1 como de PC2 en un
# grid 2x1 (patchwork)
plot_pc1 + plot_pc2 + plot_layout(ncol = 2, nrow = 1)
GRÁFICAS PCA
Utilizamos el paquete factoextra para realizar las
gráficas de los 2 primeros componentes, en concreto cla función
fviz_pca_ind, que además nos permite añadir las elipses de
confianza.
Entre los diferentes grupos de estudio, las elipses de confianza están centradas y solapadas lo que sugiere que, en el espacio de los primeros dos componentes principales, no hay una separación clara entre los grupos.
No es descartable la necesidad de examinar componentes adicionales o utilizar otras técnicas de análisis como métodos de Clustering.
# Utilizamos "factoextra" para hacer la gráfica de la PCA
# Hacemos la gráfica de los dos primeros componentes
fviz_pca_ind(pca_result,
col.ind = df_dataClin$study, # Colores de los puntos basados en los estudios
label = "none",
addEllipses = TRUE, # Añadir elipses de confianza
ellipse.type = "confidence", # Tipo de elipse, puede ser 'confidence' o 't'
legend.title = "Study",
palette = "jco", # Paleta de colores, 'jco' es solo una opción
title = "PCA Estudios")
Cuando representamos la PCA por inestabilidad de microsatélites vemos
que las dos elipses, las de MSI y las de MSS, están claramente separadas
lo que sugiere que ambos grupos tienen perfiles de expresión genética
distintos en las dimensiones capturadas por los componentes principales
(CP).
Como ya sabemos, estas diferencias son biológicamente relevantes; los
componentes que muestran esta separación probablemente están capturando
las variaciones genéticas ya conocidas de antemano.
# Hacemos la gráfica de los dos primeros componentes en relación a los grupos
fviz_pca_ind(pca_result,
col.ind = df_dataClin$msi_imputed, # Colores de los puntos basados en MSI/MSS
label = "none",
addEllipses = TRUE, # Añadir elipses de confianza
ellipse.type = "confidence", # Tipo de elipse, puede ser 'confidence' o 't'
legend.title = "MSI/MSS",
palette = "jco", # Paleta de colores, 'jco' es solo una opción
title = "PCA MSI/MSS")
Entre los dos grupos de estadio, II y III, las elipses de confianza están más o menos centradas y algo solapadas lo que sugiere que, no hay una separación clara entre los grupos. Algo un tanto previsible ya que el estadio del cáncer depende del tiempo de progreso de la enfermedad.
# Hacemos la gráfica de los dos primeros componentes en relación a los grupos
fviz_pca_ind(pca_result,
col.ind = df_dataClin$stage, # Colores de los puntos basados en Estadio
label = "none",
addEllipses = TRUE, # Añadir elipses de confianza
ellipse.type = "confidence", # Tipo de elipse, puede ser 'confidence' o 't'
legend.title = "Estadio",
palette = "jco", # Paleta de colores, 'jco' es solo una opción
title = "PCA Estadio II/III")
Se podría decir que esta gráfica és una combinación de las dos anteriores; fijándonos en las elipses de confianza, se observan dos grupos separados claramente de otros dos, debido principalmente a la contribución de la instabilidad de microsatélites y no tanto a la del estadio del CRC.
# Hacemos la gráfica de los dos primeros componentes en relación a los grupos
fviz_pca_ind(pca_result,
col.ind = df_dataClin$grupo, # Colores de los puntos basados en Grupos
label = "none",
addEllipses = TRUE, # Añadir elipses de confianza
ellipse.type = "confidence", # Tipo de elipse, puede ser 'confidence' o 't'
legend.title = "Grupo",
palette = "jco", # Paleta de colores, 'jco' es solo una opción
title = "PCA Grupos")
Al igual que con la inestabilidad de microsatélites, cuando
representamos la PCA por clasificación CMS vemos que las cuatro elipses
están claramente separadas; los 4 grupos tienen perfiles de expresión
genética distintos en las dimensiones capturadas por los CP.
Estas diferencias son biológicamente relevantes dado que los componentes
que muestran esta separación están capturando variaciones genéticas ya
conocidas.
# Hacemos la gráfica de los dos primeros componentes en relación a los grupos
fviz_pca_ind(pca_result,
col.ind = df_dataClin$cms, # Colores de los puntos basados en CMS
label = "none",
addEllipses = TRUE, # Añadir elipses de confianza
ellipse.type = "confidence", # Tipo de elipse, puede ser 'confidence' o 't'
legend.title = "CMS",
palette = "jco", # Paleta de colores, 'jco' es solo una opción
title = "PCA por CMS")
# Hacemos la gráfica de los dos primeros componentes en relación a los grupos
fviz_pca_ind(pca_result,
col.ind = df_dataClin$cms_stage, # Colores de los puntos basados en CMS
label = "none",
addEllipses = TRUE, # Añadir elipses de confianza
ellipse.type = "confidence", # Tipo de elipse, puede ser 'confidence' o 't'
legend.title = "cms_stage",
palette = "jco", # Paleta de colores, 'jco' es solo una opción
title = "PCA por CMS-Estadio")
Analizamos los datos viendo la infiltración celular.
Utilizamos el paquete MCPcounter (Becht et al. 2016) que estima la abundancia
poblacional de poblaciones de células inmunes (8 tipos) y estromales (2
tipos) que se infiltran en tejidos mediante la expresión genética.
En concreto, utilizatemos MCPcounter para estimar la
infiltración celular del estroma y el sistema inmunitario en el
microambiente tumoral (TME).
# Cargamos la librería 'MCPcounter'
require(MCPcounter)
# Con todos los datos de la matriz de expresión 26/03/2024
exmat_MCP <- MCPcounter.estimate(expression = df_exmat, featuresType = "HUGO_symbols")
# Transponemos
exmat_MCP <- t_df(exmat_MCP)
# Añadimos los datos clínicos a los datos de la matriz generada con MCPcounter
# 26/03/2024
exmat_MCP_Clin <- cbind(df_dataClin, exmat_MCP)
# Guardamos 'df_dataClin', 'df_exmat' y 'exmat_MCP_Clin' en un archivo .Rdata.
# Lo hago para poder utilizar estos datos ya calculados en otros .rmd
# 21/04/2024 save(df_dataClin, df_exmat, exmat_MCP_Clin, file =
# 'tfm_SKF.RData')
A modo de recordatorio didáctico de Inferencia Estadística:
| \(H_0\) es verdadera | \(H_0\) es falsa | |
|---|---|---|
| Rechazar \(H_0\) | Error Tipo I (α) | Decisión Correcta |
| No rechazar \(H_0\) | Decisión Correcta | Error Tipo II (β) |
El p-valor es la probabilidad del resultado del estadístico de contraste observado o de uno más alejado cuando la hipótesis nula es cierta.
El p-valor asociado a una observación del estadístico de contraste es el menor nivel de significación que nos permite rechazar la hipótesis nula.
Si el p-valor es inferior al nivel de significación \(\alpha\), rechazaremos la hipótesis nula, en caso contrario la aceptaremos.
En nuestro caso utilizamos un nivel de significación \(\alpha = 0.05\)
Para hacer una comparativa entre los 4 grupos, utilizaremos el test no paramétrico de Kruskal-Wallis. Se suele utilizar para determinar si hay diferencias significativas entre tres o más grupos independientes cuando los datos no cumplen los supuestos necesarios para realizar un ANOVA: normalidad y homocedasticidad (homogeneidad de varianzas).
# Calculamos los test no paraméticos Kruskal-Wallis (asumimos que no tenemos
# normalidad) 16/04/2024
# tests de Kruskal-Wallis para NK
kw_NK_result <- kruskal.test(exmat_MCP_Clin$`NK cells` ~ grupo, data = exmat_MCP_Clin)
# Para el resto de tipos celulares utilizamos una lista con los nombres
tipos_celulares <- c("NK cells", "T cells", "CD8 T cells", "Cytotoxic lymphocytes",
"B lineage", "Monocytic lineage", "Myeloid dendritic cells", "Neutrophils", "Endothelial cells",
"Fibroblasts")
# Aplicamos el test de Kruskal-Wallis a cada tipo celular con 'lapply'
resultados_kw <- lapply(tipos_celulares, function(cell_type) {
kruskal.test(reformulate("grupo", response = cell_type), data = exmat_MCP_Clin)
})
# Añadimos los nombres de los tipos celulares
names(resultados_kw) <- tipos_celulares
# Extraemos los valores de p de todos los tests realizados
p_valores <- sapply(resultados_kw, function(x) x$p.value)
La hipótesis nula del test de Kruskal-Wallis \(H_0\) establece que todas las poblaciones (grupos) tienen distribuciones idénticas, es decir, las medianas de todos los grupos son iguales.
\[ H_0: \forall i, j \, (i \neq j) \Rightarrow \mu_i = \mu_j \]
La hipótesis alternativa \(H_1\) establece que al menos uno de los grupos tiene una distribución diferente en cuanto a la mediana, comparada con las otras poblaciones sin especificar cuál de los grupos es diferente.
\[ H_1: \exists i \neq j \mid \mu_i \neq \mu_j \]
# Generamos todos los boxplot 27/03/2024
grupo <- exmat_MCP_Clin$grupo
# NK
bp1 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`NK cells`, grupo,
"NK cells", p_valores["NK cells"])
# T cells
bp2 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`T cells`, grupo,
"T cells", p_valores["T cells"])
# CD8 T cells
bp3 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`CD8 T cells`, grupo,
"CD8 T cells", p_valores["CD8 T cells"])
# Cytotoxic lymphocytes
bp4 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`Cytotoxic lymphocytes`,
grupo, "Cytotoxic lymphocytes", p_valores["Cytotoxic lymphocytes"])
# B lineage
bp5 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`B lineage`, grupo,
"B lineage", p_valores["B lineage"])
# Monocytic lineage
bp6 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`Monocytic lineage`,
grupo, "Monocytic lineage", p_valores["Monocytic lineage"])
# Myeloid dendritic cells
bp7 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`Myeloid dendritic cells`,
grupo, "Myeloid dendritic cells", p_valores["Myeloid dendritic cells"])
# Neutrophils
bp8 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$Neutrophils, grupo,
"Neutrophils", p_valores["Neutrophils"])
# Endothelial cells
bp9 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`Endothelial cells`,
grupo, "Endothelial cells", p_valores["Endothelial cells"])
# Fibroblasts
bp10 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$Fibroblasts, grupo,
"Fibroblasts", p_valores["Fibroblasts"])
# En un grid 3x2 (patchwork)
bp1 + bp2 + bp3 + plot_layout(ncol = 3, nrow = 1)
bp4 + bp5 + bp6 + plot_layout(ncol = 3, nrow = 1)
bp7 + bp8 + bp9 + plot_layout(ncol = 3, nrow = 1)
bp10 + plot_layout(ncol = 3, nrow = 1)
# Todos juntos para el informe del TFM 11/06/2024
bp1 + bp2 + bp3 + bp4 + bp5 + bp6 + bp7 + bp8 + bp9 + bp10 + plot_layout(ncol = 4,
nrow = 3)
NORMALIDAD Y HOMOCEDASTICIDAD
Primero compromabos la normalidad de los grupos para las NK cells.
Tabla con los resultados de los tests de Shapiro
para normalidad.
Si p-value < 0.05 rechazamos la hipótesis nula de igualdad de
estadísticos.
| Shapiro test | MSI | MSS | II | III |
|---|---|---|---|---|
| p-valor | 3.339e-05 | 1.058e-10 | 1.371e-10 | 1.089e-04 |
Dado que no hay normalidad no hace falta comprobar homocedasticidad.
Para realizar una comparativa entre dos grupos en los que asumimos no normalidad , en nuestro caso MSI vs MSS y estadio II vs estadio III, utilizamos el test no paramétrico de Mann-Whitney U también conocido como el test de Wilcoxon en el que se asume que los datos no necesitan ser normalmente distribuidos aunque deben poder ordenarse.
# Para dos grupos utilizamos el test no paramético Mann-Whitney U (asumimos que
# no tenemos normalidad) 16/04/2024
# Realizamos el test de Mann-Whitney U para los grupos MSI y MSS
mwU_NK_msi_test <- wilcox.test(grupo_msi, grupo_mss, alternative = "two.sided")
# Realizamos el test de Mann-Whitney U para los estadios II y III
mwU_NK_stage_test <- wilcox.test(grupo_sII, grupo_sIII, alternative = "two.sided")
La hipótesis nula \(H_0\) para el test de Mann-Whitney U establece que no hay diferencia en las medianas (o, más generalmente, en las distribuciones) de los dos grupos: \[ H_0: \mu_{grupo1} = \mu_{grupo2} \] La hipótesis alternativa \(H_1\) establece que hay una diferencia entre las distribuciones de los dos grupos. En nuestro caso la planteamos Bilateral (Two-sided):
\[ H_1: \mu_{grupo1} \neq \mu_{grupo2} \]
# Generamos todos los boxplot referentes a NKs 27/03/2024 Diagrama de Barras de
# NKs MSI vs MSS
bp11 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`NK cells`,
exmat_MCP_Clin$msi_imputed, "NK cells MSI vs MSS", mwU_NK_msi_test$p.value)
# Diagrama de Barras de NKs II vs III
bp12 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`NK cells`,
exmat_MCP_Clin$stage, "NK cells II vs III", mwU_NK_stage_test$p.value)
# En un grid 2x2 (patchwork)
p4 + bp1 + p1 + bp11 + p3 + bp12 + plot_layout(ncol = 2, nrow = 3)
# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp13 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`T cells`,
exmat_MCP_Clin$msi_imputed, "T cells MSI vs MSS")
# Diagrama de Barras II vs III
bp14 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`T cells`,
exmat_MCP_Clin$stage, "T cells II vs III")
# En un grid 1x3 (patchwork)
bp2 + bp13 + bp14 + plot_layout(ncol = 3, nrow = 1)
# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp15 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`CD8 T cells`,
exmat_MCP_Clin$msi_imputed, "CD8 T cells MSI vs MSS")
# Diagrama de Barras II vs III
bp16 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`CD8 T cells`,
exmat_MCP_Clin$stage, "CD8 T cells II vs III")
# En un grid 1x3 (patchwork)
bp3 + bp15 + bp16 + plot_layout(ncol = 3, nrow = 1)
# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp17 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`Cytotoxic lymphocytes`,
exmat_MCP_Clin$msi_imputed, "Cytotoxic lymphocytes MSI vs MSS")
# Diagrama de Barras II vs III
bp18 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`Cytotoxic lymphocytes`,
exmat_MCP_Clin$stage, "Cytotoxic lymphocytes II vs III")
# En un grid 1x3 (patchwork)
bp4 + bp17 + bp18 + plot_layout(ncol = 3, nrow = 1)
# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp19 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`B lineage`,
exmat_MCP_Clin$msi_imputed, "B lineage MSI vs MSS")
# Diagrama de Barras II vs III
bp20 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`B lineage`,
exmat_MCP_Clin$stage, "B lineage II vs III")
# En un grid 1x3 (patchwork)
bp5 + bp19 + bp20 + plot_layout(ncol = 3, nrow = 1)
# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp21 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`Monocytic lineage`,
exmat_MCP_Clin$msi_imputed, "Monocytic lineage MSI vs MSS")
# Diagrama de Barras II vs III
bp22 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`Monocytic lineage`,
exmat_MCP_Clin$stage, "Monocytic lineage II vs III")
# En un grid 1x3 (patchwork)
bp6 + bp21 + bp22 + plot_layout(ncol = 3, nrow = 1)
# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp23 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`Myeloid dendritic cells`,
exmat_MCP_Clin$msi_imputed, "Myeloid dendritic cells MSI vs MSS")
# Diagrama de Barras II vs III
bp24 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`Myeloid dendritic cells`,
exmat_MCP_Clin$stage, "Myeloid dendritic cells II vs III")
# En un grid 1x3 (patchwork)
bp7 + bp23 + bp24 + plot_layout(ncol = 3, nrow = 1)
# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp25 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$Neutrophils,
exmat_MCP_Clin$msi_imputed, "Neutrophils MSI vs MSS")
# Diagrama de Barras II vs III
bp26 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$Neutrophils,
exmat_MCP_Clin$stage, "Neutrophils II vs III")
# En un grid 1x3 (patchwork)
bp8 + bp25 + bp26 + plot_layout(ncol = 3, nrow = 1)
# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp27 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`Endothelial cells`,
exmat_MCP_Clin$msi_imputed, "Endothelial cells MSI vs MSS")
# Diagrama de Barras II vs III
bp28 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`Endothelial cells`,
exmat_MCP_Clin$stage, "Endothelial cells II vs III")
# En un grid 1x3 (patchwork)
bp9 + bp27 + bp28 + plot_layout(ncol = 3, nrow = 1)
# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp29 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$Fibroblasts,
exmat_MCP_Clin$msi_imputed, "Fibroblasts MSI vs MSS")
# Diagrama de Barras II vs III
bp30 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$Fibroblasts,
exmat_MCP_Clin$stage, "Fibroblasts II vs III")
# En un grid 1x3 (patchwork)
bp10 + bp29 + bp30 + plot_layout(ncol = 3, nrow = 1)
# En un grid 1x3
bp1 + bp11 + bp12 + plot_layout(ncol = 3, nrow = 1)
La clasificación CMS divide los tumores de CRC en cuatro subtipos principales, basados en características moleculares distintas, lo que refleja diferencias en la biología del tumor, la interacción con el TME, el pronóstico y la respuesta a tratamientos específicos. (Roelands et al. 2017)
| Subtipo | Características | Inmunidad | Prevalencia |
|---|---|---|---|
| CMS1 (Subtipo Inmuno) | Alta inmunogenicidad, frecuentemente asociada con MSI y una alta carga mutacional. | Tienen una infiltración activa de células inmunitarias efectoras, lo que puede hacerlos más susceptibles a terapias inmunogénicas. | Aproximadamente el 14% de los casos de CRC. |
| CMS2 (Subtipo Canónico) | Tumores con alta expresión de genes relacionados con el ciclo celular y la señalización WNT/MYC, presentando una marcada proliferación celular. | Son generalmente poco inmunogénicos, con menor infiltración inmune comparado con CMS1 y CMS4. | Es el subtipo más común, abarcando alrededor del 37% de los casos de CRC. |
| CMS3 (Subtipo Metabólico) | Presenta alteraciones metabólicas, incluyendo un metabolismo energético distintivo, con una menor carga mutacional en comparación con CMS1. | Similar a CMS2, muestra baja inmunogenicidad. | Constituye aproximadamente el 13% de los casos de CRC |
| CMS4 (Subtipo Mesenquimal) | Se caracteriza por la activación de vías relacionadas con la angiogénesis, la activación de células estromales y la señalización TGF-β, lo que conduce a un entorno inmunosupresor. | Aunque presenta una significativa infiltración inmune, esta es predominantemente de naturaleza supresora, con una alta presencia de células estromales. | Aproximadamente el 23% de los casos de CRC caen en este subtipo. |
| Vemos la infiltración celular según CMS e inestabilidad de microsatélites.
# Calculamos los test no paraméticos Kruskal-Wallis (asumimos que no tenemos
# normalidad) 16/04/2024
# Aplicamos el test de Kruskal-Wallis a cada tipo celular con 'lapply'
resultados_kw_cms <- lapply(tipos_celulares, function(cell_type) {
kruskal.test(reformulate("cms", response = cell_type), data = exmat_MCP_Clin)
})
# Añadimos los nombres de los tipos celulares
names(resultados_kw_cms) <- tipos_celulares
# Extraemos los valores de p de todos los tests realizados
p_valores_cms <- sapply(resultados_kw_cms, function(x) x$p.value)
# Quitamos los NA's de la columna cms
exmat_MCP_Clin_CMS <- exmat_MCP_Clin[!is.na(exmat_MCP_Clin$cms), ]
# Generamos todos los boxplot 02/04/2024
cms <- exmat_MCP_Clin_CMS$cms
# NK
bp1 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`NK cells`,
cms, "NK cells", p_valores_cms["NK cells"])
# T cells
bp2 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`T cells`,
cms, "T cells", p_valores_cms["T cells"])
# CD8 T cells
bp3 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`CD8 T cells`,
cms, "CD8 T cells", p_valores_cms["CD8 T cells"])
# Cytotoxic lymphocytes
bp4 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`Cytotoxic lymphocytes`,
cms, "Cytotoxic lymphocytes", p_valores_cms["Cytotoxic lymphocytes"])
# B lineage
bp5 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`B lineage`,
cms, "B lineage", p_valores_cms["B lineage"])
# Monocytic lineage
bp6 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`Monocytic lineage`,
cms, "Monocytic lineage", p_valores_cms["Monocytic lineage"])
# Myeloid dendritic cells
bp7 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`Myeloid dendritic cells`,
cms, "Myeloid dendritic cells", p_valores_cms["Myeloid dendritic cells"])
# Neutrophils
bp8 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$Neutrophils,
cms, "Neutrophils", p_valores_cms["Neutrophils"])
# Endothelial cells
bp9 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`Endothelial cells`,
cms, "Endothelial cells", p_valores_cms["Endothelial cells"])
# Fibroblasts
bp10 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$Fibroblasts,
cms, "Fibroblasts", p_valores_cms["Fibroblasts"])
# En un grid 3x2 (patchwork)
p2 + bp1 + bp2 + plot_layout(ncol = 3, nrow = 1)
bp3 + bp4 + bp5 + plot_layout(ncol = 3, nrow = 1)
bp6 + bp7 + bp8 + plot_layout(ncol = 3, nrow = 1)
bp9 + bp10 + plot_layout(ncol = 3, nrow = 1)
Los tumores CMS2 y CMS3 son poco inmunogénicos; CMS1 y CMS4 difieren
en el tipo de infiltración inmune.
Los tumores CMS1 tienden a acumular una gran cantidad de mutaciones
debido al estado de MSI y atraen células efectoras inmunes.
Los tumores CMS4 exhiben un TME inmunosupresor con infiltración
estromal. (Lanuza et al. 2022)
| Vemos la infiltración celular según CMS y estadio.
# Calculamos los test no paraméticos Kruskal-Wallis (asumimos que no tenemos
# normalidad) 22/04/2024
# Aplicamos el test de Kruskal-Wallis a cada tipo celular con 'lapply'
resultados_kw_cms_stage <- lapply(tipos_celulares, function(cell_type) {
kruskal.test(reformulate("cms_stage", response = cell_type), data = exmat_MCP_Clin)
})
# Añadimos los nombres de los tipos celulares
names(resultados_kw_cms_stage) <- tipos_celulares
# Extraemos los valores de p de todos los tests realizados
p_valores_cms_stage <- sapply(resultados_kw_cms_stage, function(x) x$p.value)
# Quitamos los NA's de la columna cms
exmat_MCP_Clin_CMS_stage <- exmat_MCP_Clin[!is.na(exmat_MCP_Clin$cms_stage), ]
# Generamos todos los boxplot 02/04/2024
cms_stage <- exmat_MCP_Clin_CMS$cms_stage
# NK
bp1 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS_stage, cms_stage, exmat_MCP_Clin_CMS_stage$`NK cells`,
cms_stage, "NK cells", p_valores_cms_stage["NK cells"])
# T cells
bp2 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS_stage, cms_stage, exmat_MCP_Clin_CMS_stage$`T cells`,
cms_stage, "T cells", p_valores_cms_stage["T cells"])
# CD8 T cells
bp3 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms_stage, exmat_MCP_Clin_CMS$`CD8 T cells`,
cms_stage, "CD8 T cells", p_valores_cms_stage["CD8 T cells"])
# Cytotoxic lymphocytes
bp4 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS_stage, cms_stage, exmat_MCP_Clin_CMS_stage$`Cytotoxic lymphocytes`,
cms_stage, "Cytotoxic lymphocytes", p_valores_cms_stage["Cytotoxic lymphocytes"])
# B lineage
bp5 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS_stage, cms_stage, exmat_MCP_Clin_CMS_stage$`B lineage`,
cms_stage, "B lineage", p_valores_cms_stage["B lineage"])
# Monocytic lineage
bp6 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS_stage, cms_stage, exmat_MCP_Clin_CMS_stage$`Monocytic lineage`,
cms_stage, "Monocytic lineage", p_valores_cms_stage["Monocytic lineage"])
# Myeloid dendritic cells
bp7 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS_stage, cms_stage, exmat_MCP_Clin_CMS_stage$`Myeloid dendritic cells`,
cms_stage, "Myeloid dendritic cells", p_valores_cms_stage["Myeloid dendritic cells"])
# Neutrophils
bp8 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS_stage, cms_stage, exmat_MCP_Clin_CMS_stage$Neutrophils,
cms_stage, "Neutrophils", p_valores_cms_stage["Neutrophils"])
# Endothelial cells
bp9 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS_stage, cms_stage, exmat_MCP_Clin_CMS_stage$`Endothelial cells`,
cms_stage, "Endothelial cells", p_valores_cms_stage["Endothelial cells"])
# Fibroblasts
bp10 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS_stage, cms_stage, exmat_MCP_Clin_CMS_stage$Fibroblasts,
cms_stage, "Fibroblasts", p_valores_cms_stage["Fibroblasts"])
# En un grid 2x1 (patchwork)
bp1 + bp2 + plot_layout(ncol = 2, nrow = 1)
bp3 + bp4 + plot_layout(ncol = 2, nrow = 1)
bp5 + bp6 + plot_layout(ncol = 2, nrow = 1)
bp7 + bp8 + plot_layout(ncol = 2, nrow = 1)
bp9 + bp10 + plot_layout(ncol = 2, nrow = 1)
p5 + plot_layout(ncol = 2, nrow = 1)
Hacemos un data frame con los datos de expresión de los genes del Complejo Mayor de Hitocompatibilidad I HLA-A, HLA-B, HLA-C y el gen marcador de las células NK, NKp46 (NCR1 en nomenclatura SYMBOL). Además añadimos la columna HLA-ABC con la suma de expresión de los 3 genes.
# Nos quedamos solo con los genes del complejo mayor de hitocompatibilidad I
# HLA-ABC y con el marcador de células NK, el gen NKp46 (NCR1 en nomenclatura
# SYMBOL)
hla_nk_genes <- c("HLA-A", "HLA-B", "HLA-C", "HLA-E", "NCR1")
df_hla_nk_expression <- df_exmat_t[hla_nk_genes]
# Calculamos la mediana de la columna HLA-A
median_hla_a <- median(df_hla_nk_expression$`HLA-A`, na.rm = TRUE)
# Calculamos la mediana de la columna HLA-A
median_hla_b <- median(df_hla_nk_expression$`HLA-B`, na.rm = TRUE)
# Calculamos la mediana de la columna HLA-A
median_hla_c <- median(df_hla_nk_expression$`HLA-C`, na.rm = TRUE)
# Calculamos la mediana de la columna HLA-E
median_hla_e <- median(df_hla_nk_expression$`HLA-E`, na.rm = TRUE)
# Añadimos una nueva columna que indique si el nivel de HLA-A es 'Low' o 'High'
df_hla_nk_expression$HLA_A_cat <- ifelse(df_hla_nk_expression$`HLA-A` < median_hla_a,
"Low", "High")
# Convertimos la columna 'HLA_A_cat' a factor
df_hla_nk_expression$HLA_A_cat <- factor(df_hla_nk_expression$HLA_A_cat, levels = c("Low",
"High"))
# Añadimos una nueva columna que indique si el nivel de HLA-B es 'Low' o 'High'
df_hla_nk_expression$HLA_B_cat <- ifelse(df_hla_nk_expression$`HLA-B` < median_hla_b,
"Low", "High")
# Convertimos la columna 'HLA_B_cat' a factor
df_hla_nk_expression$HLA_B_cat <- factor(df_hla_nk_expression$HLA_B_cat, levels = c("Low",
"High"))
# Añadimos una nueva columna que indique si el nivel de HLA-C es 'Low' o 'High'
df_hla_nk_expression$HLA_C_cat <- ifelse(df_hla_nk_expression$`HLA-C` < median_hla_c,
"Low", "High")
# Convertimos la columna 'HLA_C_cat' a factor
df_hla_nk_expression$HLA_C_cat <- factor(df_hla_nk_expression$HLA_C_cat, levels = c("Low",
"High"))
# Añadimos una nueva columna que indique si el nivel de HLA-E es 'Low' o 'High'
df_hla_nk_expression$HLA_E_cat <- ifelse(df_hla_nk_expression$`HLA-E` < median_hla_e,
"Low", "High")
# Convertimos la columna 'HLA_E_cat' a factor
df_hla_nk_expression$HLA_E_cat <- factor(df_hla_nk_expression$HLA_E_cat, levels = c("Low",
"High"))
# Añadimos los datos clínicos
df_hla_nk_expression <- cbind(exmat_MCP_Clin, df_hla_nk_expression)
Boxplot comparando la expresión del gen NKp46 entre las categorías
baja (Low) y alta (High) de HLA-A. Para realizar las categorías se ha
utilizado la mediana de la expresión génica como punto de corte.
Se ha estratificando por estadio y MSI/MSS. Utilizamos el test no
paramétrico de Wilcoxon para evaluar las diferencias.
# Función para generar gráficos de caja y pruebas de Mann-Whitney U
boxplot_MWU <- function(data, stage, msi, hla_var, nk_var, title_prefix) {
# Filtrar los datos
df_filtered <- sqldf(paste0("SELECT * FROM df_hla_nk_expression WHERE stage = '",
stage, "' AND msi_imputed = '", msi, "'"))
# Realizar el test de Mann-Whitney U para los grupos Low y High
nLow <- df_filtered$NCR1[df_filtered[[hla_var]] == "Low"]
nHigh <- df_filtered$NCR1[df_filtered[[hla_var]] == "High"]
mwU_test <- wilcox.test(nLow, nHigh, alternative = "two.sided")
# Generar el título completo
title <- paste0(title_prefix, ", ", stage, ", ", msi)
# Generar el boxplot
bp <- boxplot_MCPcounter(data = df_filtered, x = df_filtered[[hla_var]], y = df_filtered[[nk_var]],
fill = df_filtered[[hla_var]], title = title, p_value = mwU_test$p.value,
xlabel = hla_var, ylabel = nk_var)
return(bp)
}
# Generar los gráficos para cada combinación de stage, msi y HLA variable
bp35 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSI", hla_var = "HLA_A_cat",
nk_var = "NCR1", title_prefix = "StageII")
bp36 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSS", hla_var = "HLA_A_cat",
nk_var = "NCR1", title_prefix = "StageII")
bp37 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSI", hla_var = "HLA_A_cat",
nk_var = "NCR1", title_prefix = "StageIII")
bp38 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSS", hla_var = "HLA_A_cat",
nk_var = "NCR1", title_prefix = "StageIII")
# Visualizar los gráficos en un grid 2x2 (patchwork)
bp35 + bp36 + bp37 + bp38 + plot_layout(ncol = 2, nrow = 2)
No da significación en ningún caso. No hay una diferencia signficativa en la mediana de la expresión del gen NCR1 entre los grupos con baja y alta expresión de HLA-A.
# Seleccionamos estadio tipo II
df_hla_nk_expression_II <- sqldf('SELECT * FROM df_hla_nk_expression WHERE stage = "II"')
# Creamos el gráfico de dispersión para Estadio II
pd2 <- ggplot(df_hla_nk_expression_II, aes(x = `HLA-A`, y = NCR1)) +
geom_point(aes(color = msi_imputed, shape = msi_imputed), size = 3, stroke = 1) +
labs(title = "Estadio II",
x = "HLA-A",
y = "NKP46 (NCR1)",
color = "MSI/MSS",
shape = "MSI/MSS") +
scale_shape_manual(values = c(16, 2)) + # cuadrado lleno, triángulo
scale_color_manual(values = c("darkorange2", "dodgerblue4")) +
theme_classic() +
theme(
axis.line = element_line(color = "black"),
legend.position = c(0.1, 0.9))
# Añadimos líneas de regresión lineal para cada grupo MSI y MSS
pd2 <- pd2 + geom_smooth(method = "lm", aes(color = msi_imputed), se = FALSE)
# Visualizamos el gráfico
pd2
# Ajustamos un modelo lineal para Estadio II
lmod_ser <- lm(NCR1 ~ `HLA-A` * msi_imputed, data = df_hla_nk_expression_II)
# Hacemos una contraste de paralelismo de rectas con anova
anova(lmod_ser)
Analysis of Variance Table
Response: NCR1
Df Sum Sq Mean Sq F value Pr(>F)
`HLA-A` 1 0.16 0.1625 0.2164 0.641981
msi_imputed 1 1.49 1.4871 1.9799 0.159859
`HLA-A`:msi_imputed 1 6.93 6.9312 9.2282 0.002474 **
Residuals 680 510.74 0.7511
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
El efecto principal de HLA-A sobre NCR1 no es significativo (p = 0.641981). Esto sugiere que no hay suficiente evidencia para concluir que HLA-A tiene un efecto significativo en NCR1 por sí solo.
El efecto msi_imputed sobre NCR1 no es significativo (p = 0.159859). Esto indica que no hay suficiente evidencia para concluir que msi_imputed tiene un efecto significativo en NCR1 por sí solo.
Vemos que hay significación en la hipótesis de paralelismo (Pr: 0.002474).La interacción entre HLA-A y msi_imputed es significativa. Esto indica que el efecto de HLA-A sobre NCR1 depende de los valores de msi_imputed. Dado que el p-valor es menor que 0.01 (**), hay una fuerte evidencia contra la hipótesis nula de que no hay interacción entre HLA-A y msi_imputed.
Aquí se observa para los tumores MSI una relación inversa entre la expresión de NCR1 y HLA-A; Cuanta menor expresión de HLA’s hay una mayor expresión de NKp46, es decir, cuanta menor expresión de HLA-A hay una mayor infiltración de células NK activadas.
# Seleccionamos estadio tipo II
df_hla_nk_expression_III <- sqldf('SELECT * FROM df_hla_nk_expression WHERE stage = "III"')
# Creamos el gráfico de dispersión para Estadio II
pd3 <- ggplot(df_hla_nk_expression_III, aes(x = `HLA-A`, y = NCR1)) +
geom_point(aes(color = msi_imputed, shape = msi_imputed), size = 3, stroke = 1) +
labs(title = "Estadio III",
x = "HLA-A",
y = "NKP46 (NCR1)",
color = "MSI/MSS",
shape = "MSI/MSS") +
scale_shape_manual(values = c(16, 2)) + # cuadrado lleno, triángulo
scale_color_manual(values = c("darkorange2", "dodgerblue4")) +
theme_classic() +
theme(
axis.line = element_line(color = "black"),
legend.position = c(0.1, 0.9))
# Añadimos líneas de regresión lineal para cada grupo MSI y MSS
pd3 <- pd3 + geom_smooth(method = "lm", aes(color = msi_imputed), se = FALSE)
# Visualizamos el gráfico
pd3
# Ajustamos un modelo lineal para Estadio II
lmod_ser <- lm(NCR1 ~ `HLA-A` * msi_imputed, data = df_hla_nk_expression_III)
# Hacemos una contraste de paralelismo de rectas con anova
anova(lmod_ser)
Analysis of Variance Table
Response: NCR1
Df Sum Sq Mean Sq F value Pr(>F)
`HLA-A` 1 3.172 3.1717 4.8625 0.028028 *
msi_imputed 1 5.148 5.1481 7.8925 0.005213 **
`HLA-A`:msi_imputed 1 0.036 0.0356 0.0545 0.815465
Residuals 391 255.037 0.6523
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
El efecto principal de HLA-A sobre NCR1 es significativo (p = 0.028028). Esto sugiere que hay una relación significativa entre HLA-A y NCR1 por sí solo. Dado que el p-valor es menor que 0.05 (*), hay suficiente evidencia para rechazar la hipótesis nula de que HLA-A no tiene un efecto sobre NCR1.
El efecto principal de msi_imputed sobre NCR1 es significativo (p = 0.005213). Esto indica que hay una relación significativa entre msi_imputed y NCR1 por sí solo. Dado que el p-valor es menor que 0.01 (**), hay una fuerte evidencia para rechazar la hipótesis nula de que msi_imputed no tiene un efecto sobre NCR1.
La interacción entre HLA-A y msi_imputed no es significativa (p = 0.815465). Esto indica que no hay suficiente evidencia para concluir que el efecto de HLA-A sobre NCR1 depende de los valores de msi_imputed.
En este caso no hay significación en la hipótesis de paralelismo (Pr: 0.738591) aunque sí en la de coincidencia de las rectas (Pr: 0.006508). Resultado también de difícil interpretación: los efectos principales de HLA-ABC y msi_imputed son ambos estadísticamente significativos y aportan individualmente a la explicación de la variabilidad en NCR1, con msi_imputed mostrando un efecto algo más fuerte según valor p más bajos.
Aquí se intuye una relación inversa entre la expresión de NCR1 y HLA-A; es decir, una relación inversa entre infiltración de células NK y expresión de HLA’s.
Boxplot comparando la expresión del gen NKp46 entre las categorías
baja (Low) y alta (High) de HLA-B. Para realizar las categorías se ha
utilizado la mediana de la expresión génica como punto de corte.
Se ha estratificando por estadio y MSI/MSS. Utilizamos el test no
paramétrico de Wilcoxon para evaluar las diferencias.
# Generar los gráficos para cada combinación de stage, msi y HLA variable
bp39 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSI", hla_var = "HLA_B_cat",
nk_var = "NCR1", title_prefix = "StageII")
bp40 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSS", hla_var = "HLA_B_cat",
nk_var = "NCR1", title_prefix = "StageII")
bp41 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSI", hla_var = "HLA_B_cat",
nk_var = "NCR1", title_prefix = "StageIII")
bp42 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSS", hla_var = "HLA_B_cat",
nk_var = "NCR1", title_prefix = "StageIII")
# Visualizar los gráficos en un grid 2x2 (patchwork)
bp39 + bp40 + bp41 + bp42 + plot_layout(ncol = 2, nrow = 2)
Nos da significación en el primer caso StageII, MSI.
Boxplot comparando la expresión del gen NKp46 entre las categorías
baja (Low) y alta (High) de HLA-C. Para realizar las categorías se ha
utilizado la mediana de la expresión génica como punto de corte.
Se ha estratificando por estadio y MSI/MSS. Utilizamos el test no
paramétrico de Wilcoxon para evaluar las diferencias.
# Generar los gráficos para cada combinación de stage, msi y HLA variable
bp43 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSI", hla_var = "HLA_C_cat",
nk_var = "NCR1", title_prefix = "StageII")
bp44 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSS", hla_var = "HLA_C_cat",
nk_var = "NCR1", title_prefix = "StageII")
bp45 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSI", hla_var = "HLA_C_cat",
nk_var = "NCR1", title_prefix = "StageIII")
bp46 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSS", hla_var = "HLA_C_cat",
nk_var = "NCR1", title_prefix = "StageIII")
# Visualizar los gráficos en un grid 2x2 (patchwork)
bp43 + bp44 + bp45 + bp46 + plot_layout(ncol = 2, nrow = 2)
Nos da significación en los casos StageII, MSS y StageIII, MSS.
Boxplot comparando la expresión del gen NKp46 entre las categorías
baja (Low) y alta (High) de HLA-E. Para realizar las categorías se ha
utilizado la mediana de la expresión génica como punto de corte.
Se ha estratificando por estadio y MSI/MSS. Utilizamos el test no
paramétrico de Wilcoxon para evaluar las diferencias.
bp47 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSI", hla_var = "HLA_E_cat",
nk_var = "NCR1", title_prefix = "StageII")
bp48 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSS", hla_var = "HLA_E_cat",
nk_var = "NCR1", title_prefix = "StageII")
bp49 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSI", hla_var = "HLA_E_cat",
nk_var = "NCR1", title_prefix = "StageIII")
bp50 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSS", hla_var = "HLA_E_cat",
nk_var = "NCR1", title_prefix = "StageIII")
# Visualizar los gráficos en un grid 2x2
bp47 + bp48 + bp49 + bp50 + plot_layout(ncol = 2, nrow = 2)
Nos da significación en los casos StageII, MSS y StageIII, MSS.
# Seleccionamos estadio tipo II ya lo hemos hecho antes
df_hla_nk_expression_MSS <- sqldf('SELECT * FROM df_hla_nk_expression WHERE msi_imputed = "MSS"')
# Creamos el gráfico de dispersión para Estadio II
pd4 <- ggplot(df_hla_nk_expression_MSS, aes(x = `HLA-E`, y = NCR1)) +
geom_point(aes(color = stage, shape = stage), size = 3, stroke = 1) +
labs(title = "MSS",
x = "HLA-E",
y = "NKP46 (NCR1)",
color = "II/III",
shape = "II/III") +
scale_shape_manual(values = c(16, 2)) + # cuadrado lleno, triángulo
scale_color_manual(values = c("darkorange2", "dodgerblue4")) +
theme_classic() +
theme(
axis.line = element_line(color = "black"),
legend.position = c(0.1, 0.9))
# Añadimos líneas de regresión lineal para cada grupo MSI y MSS
pd4 <- pd4 + geom_smooth(method = "lm", aes(color = stage), se = FALSE)
# Visualizamos el gráfico
pd4
# Ajustamos un modelo lineal para Estadio II
lmod_ser <- lm(NCR1 ~ `HLA-A` * stage, data = df_hla_nk_expression_MSS)
# Hacemos una contraste de paralelismo de rectas con anova
anova(lmod_ser)
Analysis of Variance Table
Response: NCR1
Df Sum Sq Mean Sq F value Pr(>F)
`HLA-A` 1 5.67 5.6717 8.1903 0.004317 **
stage 1 1.23 1.2299 1.7760 0.183005
`HLA-A`:stage 1 0.00 0.0013 0.0018 0.965735
Residuals 836 578.92 0.6925
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
El efecto principal de HLA-A en NCR1 es significativo (p = 0.004317). Esto significa que hay una relación significativa entre HLA-A y la infiltración de células NK activadas. Dado que el p-valor es menor que 0.01 (**), esto indica una fuerte evidencia contra la hipótesis nula de que HLA-A no tiene un efecto sobre NCR1.
El efecto de stage en NCR1 no es significativo (p = 0.183005). Esto sugiere que no hay suficiente evidencia para concluir que el estadio tenga un efecto sobre NCR1 cuando se considera solo este factor.
Vemos que hay significación en la hipótesis de paralelismo
HLA-A:stage (Pr: 0.965735), no son paralelas. La
interacción entre HLA-A y stage no es significativa (p = 0.965735). Esto
indica que no hay suficiente evidencia para concluir que el efecto de
HLA-A sobre NCR1 depende del estadio de la muestra.
# Generar los gráficos para cada combinación de stage, msi y HLA variable
bp51 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSI", hla_var = "HLA_A_cat",
nk_var = "NK cells", title_prefix = "StageII")
bp52 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSS", hla_var = "HLA_A_cat",
nk_var = "NK cells", title_prefix = "StageII")
bp53 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSI", hla_var = "HLA_A_cat",
nk_var = "NK cells", title_prefix = "StageIII")
bp54 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSS", hla_var = "HLA_A_cat",
nk_var = "NK cells", title_prefix = "StageIII")
# Visualizar los gráficos en un grid 2x2 (patchwork)
bp51 + bp52 + bp53 + bp54 + plot_layout(ncol = 2, nrow = 2)
# Generar los gráficos para cada combinación de stage, msi y HLA variable
bp55 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSI", hla_var = "HLA_B_cat",
nk_var = "NK cells", title_prefix = "StageII")
bp56 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSS", hla_var = "HLA_B_cat",
nk_var = "NK cells", title_prefix = "StageII")
bp57 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSI", hla_var = "HLA_B_cat",
nk_var = "NK cells", title_prefix = "StageIII")
bp58 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSS", hla_var = "HLA_B_cat",
nk_var = "NK cells", title_prefix = "StageIII")
# Visualizar los gráficos en un grid 2x2 (patchwork)
bp55 + bp56 + bp57 + bp58 + plot_layout(ncol = 2, nrow = 2)
# Generar los gráficos para cada combinación de stage, msi y HLA variable
bp59 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSI", hla_var = "HLA_C_cat",
nk_var = "NK cells", title_prefix = "StageII")
bp60 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSS", hla_var = "HLA_C_cat",
nk_var = "NK cells", title_prefix = "StageII")
bp61 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSI", hla_var = "HLA_C_cat",
nk_var = "NK cells", title_prefix = "StageIII")
bp62 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSS", hla_var = "HLA_C_cat",
nk_var = "NK cells", title_prefix = "StageIII")
# Visualizar los gráficos en un grid 2x2 (patchwork)
bp59 + bp60 + bp61 + bp62 + plot_layout(ncol = 2, nrow = 2)
# Generar los gráficos para cada combinación de stage, msi y HLA variable
bp63 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSI", hla_var = "HLA_E_cat",
nk_var = "NK cells", title_prefix = "StageII")
bp64 <- boxplot_MWU(data = df_hla_nk_expression, stage = "II", msi = "MSS", hla_var = "HLA_E_cat",
nk_var = "NK cells", title_prefix = "StageII")
bp65 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSI", hla_var = "HLA_E_cat",
nk_var = "NK cells", title_prefix = "StageIII")
bp66 <- boxplot_MWU(data = df_hla_nk_expression, stage = "III", msi = "MSS", hla_var = "HLA_E_cat",
nk_var = "NK cells", title_prefix = "StageIII")
# Visualizar los gráficos en un grid 2x2 (patchwork)
bp63 + bp64 + bp65 + bp66 + plot_layout(ncol = 2, nrow = 2)
El algoritmo “Random Forest” es un método de aprendizaje supervisado
utilizado para la clasificación en el campo de machine learning (ML). Se
utilizan en una amplia gama de campos como, por ejemplo, en el
diagnóstico médico.
Aplicaremos este algoritmo a nuestros datos para clasificar los CRC así
como el estadio de los mismos.
# library(randomForest) library(caret) library(dplyr)
# Función para crear y evaluar un modelo Random Forest con 100 árboles
run_randomForest_model <- function(data, target_name = "class", p_train = 0.67, seed = 12357,
ntree = 100) {
set.seed(seed)
# Preparamos los datos
clean_data <- na.omit(data)
# Obtenemos el número total de filas (datos) n <- nrow(clean_data)
# Obtenemos los índices de filas al azar para entrenamiento indices_train
# <- sample(1:n, round(p_train * n))
# Usar createDataPartition para obtener índices de entrenamiento
indices_train <- createDataPartition(y = clean_data$class, p = params$p.train,
list = FALSE)
# Seleccionamos las filas para entrenamiento
data_train <- clean_data[indices_train, ]
# Seleccionamos las filas restantes para pruebas
data_test <- clean_data[-indices_train, ]
# Creamos los labels de class de 'train'
train_labels <- data_train[[target_name]]
# Creamos los labels de class de 'test'
test_labels <- data_test[[target_name]]
# Creamos el modelo
model <- randomForest(data_train[-1], train_labels, ntree = ntree)
# Predicción y evaluación
predictions <- predict(model, data_test[-1], type = "response")
results <- confusionMatrix(table(predictions, test_labels))
return(list(model = model, results = results))
}
# Añadimos la variable 'class' que queramos para los modelos ML
df_exmat_stage <- cbind(df_dataClin["stage"], df_exmat_t)
# Cambiamos el nombre de la primera columna a 'class'
names(df_exmat_stage)[1] <- "class"
# Ejecutamos el modelo Random Forest
model_results <- run_randomForest_model(df_exmat_stage)
# Vemos la información del modelo
model_results$model
Call:
randomForest(x = data_train[-1], y = train_labels, ntree = 100)
Type of random forest: classification
Number of trees: 100
No. of variables tried at each split: 123
OOB estimate of error rate: 37.36%
Confusion matrix:
II III class.error
II 413 43 0.09429825
III 226 38 0.85606061
# Vemos el summary del modelo summary(model_results$model) Vemos los resultados
# de la Matriz de Confusión de las predicciones (test)
print(model_results$results)
Confusion Matrix and Statistics
test_labels
predictions II III
II 215 114
III 13 17
Accuracy : 0.6462
95% CI : (0.5943, 0.6957)
No Information Rate : 0.6351
P-Value [Acc > NIR] : 0.3522
Kappa : 0.087
Mcnemar's Test P-Value : <2e-16
Sensitivity : 0.9430
Specificity : 0.1298
Pos Pred Value : 0.6535
Neg Pred Value : 0.5667
Prevalence : 0.6351
Detection Rate : 0.5989
Detection Prevalence : 0.9164
Balanced Accuracy : 0.5364
'Positive' Class : II
# Regogemos los datos para comparar de la predicción
stage_Accuracy <- model_results$results$overall["Accuracy"]
stage_Kappa <- model_results$results$overall["Kappa"]
stage_AccuracyLower <- model_results$results$overall["AccuracyLower"]
stage_AccuracyUpper <- model_results$results$overall["AccuracyUpper"]
# Añadimos la variable 'class' que queramos para los modelos ML
df_exmat_msi <- cbind(df_dataClin["msi_imputed"], df_exmat_t)
# Cambiamos el nombre de la primera columna a 'class'
names(df_exmat_msi)[1] <- "class"
# Ejecutamos el modelo Random Forest
model_results <- run_randomForest_model(df_exmat_msi, ntree = 200)
# Vemos la información del modelo
model_results$model
Call:
randomForest(x = data_train[-1], y = train_labels, ntree = ntree)
Type of random forest: classification
Number of trees: 200
No. of variables tried at each split: 123
OOB estimate of error rate: 12.64%
Confusion matrix:
MSI MSS class.error
MSI 94 66 0.41250000
MSS 25 535 0.04464286
# Vemos el summary del modelo summary(model_results$model) Vemos los resultados
# de la Matriz de Confusión de las predicciones (test)
print(model_results$results)
Confusion Matrix and Statistics
test_labels
predictions MSI MSS
MSI 34 7
MSS 45 273
Accuracy : 0.8552
95% CI : (0.8144, 0.8899)
No Information Rate : 0.7799
P-Value [Acc > NIR] : 0.0002111
Kappa : 0.49
Mcnemar's Test P-Value : 2.882e-07
Sensitivity : 0.43038
Specificity : 0.97500
Pos Pred Value : 0.82927
Neg Pred Value : 0.85849
Prevalence : 0.22006
Detection Rate : 0.09471
Detection Prevalence : 0.11421
Balanced Accuracy : 0.70269
'Positive' Class : MSI
# Regogemos los datos para comparar de la predicción
msi_Accuracy <- model_results$results$overall["Accuracy"]
msi_Kappa <- model_results$results$overall["Kappa"]
msi_AccuracyLower <- model_results$results$overall["AccuracyLower"]
msi_AccuracyUpper <- model_results$results$overall["AccuracyUpper"]
# Añadimos la variable 'class' que queramos para los modelos ML
df_exmat_msi_stage <- cbind(df_dataClin["grupo"], df_exmat_t)
# Cambiamos el nombre de la primera columna a 'class'
names(df_exmat_msi_stage)[1] <- "class"
# Ejecutamos el modelo Random Forest
model_results <- run_randomForest_model(df_exmat_msi_stage, ntree = 200)
# Vemos la información del modelo
model_results$model
Call:
randomForest(x = data_train[-1], y = train_labels, ntree = ntree)
Type of random forest: classification
Number of trees: 200
No. of variables tried at each split: 123
OOB estimate of error rate: 45.08%
Confusion matrix:
II_MSI II_MSS III_MSI III_MSS class.error
II_MSI 71 43 0 2 0.3879310
II_MSS 13 291 0 37 0.1466276
III_MSI 26 12 1 5 0.9772727
III_MSS 4 183 0 33 0.8500000
# Vemos el summary del modelo summary(model_results$model) Vemos los resultados
# de la Matriz de Confusión de las predicciones (test)
print(model_results$results)
Confusion Matrix and Statistics
test_labels
predictions II_MSI II_MSS III_MSI III_MSS
II_MSI 34 8 13 8
II_MSS 22 154 8 91
III_MSI 0 0 0 0
III_MSS 1 8 1 10
Overall Statistics
Accuracy : 0.5531
95% CI : (0.4999, 0.6053)
No Information Rate : 0.4749
P-Value [Acc > NIR] : 0.001813
Kappa : 0.2428
Mcnemar's Test P-Value : < 2.2e-16
Statistics by Class:
Class: II_MSI Class: II_MSS Class: III_MSI Class: III_MSS
Sensitivity 0.59649 0.9059 0.00000 0.09174
Specificity 0.90365 0.3564 1.00000 0.95984
Pos Pred Value 0.53968 0.5600 NaN 0.50000
Neg Pred Value 0.92203 0.8072 0.93855 0.70710
Prevalence 0.15922 0.4749 0.06145 0.30447
Detection Rate 0.09497 0.4302 0.00000 0.02793
Detection Prevalence 0.17598 0.7682 0.00000 0.05587
Balanced Accuracy 0.75007 0.6311 0.50000 0.52579
# Regogemos los datos para comparar de la predicción
msi_mss_stage_Accuracy <- model_results$results$overall["Accuracy"]
msi_mss_stage_Kappa <- model_results$results$overall["Kappa"]
msi_mss_stage_AccuracyLower <- model_results$results$overall["AccuracyLower"]
msi_mss_stage_AccuracyUpper <- model_results$results$overall["AccuracyUpper"]
# Añadimos la variable 'class' que queramos para los modelos ML
df_exmat_cms <- cbind(df_dataClin["cms"], df_exmat_t)
# Cambiamos el nombre de la primera columna a 'class'
names(df_exmat_cms)[1] <- "class"
# Ejecutamos el modelo Random Forest
model_results <- run_randomForest_model(df_exmat_cms)
# Vemos la información del modelo
model_results$model
Call:
randomForest(x = data_train[-1], y = train_labels, ntree = 100)
Type of random forest: classification
Number of trees: 100
No. of variables tried at each split: 123
OOB estimate of error rate: 11.36%
Confusion matrix:
CMS1 CMS2 CMS3 CMS4 class.error
CMS1 109 7 4 10 0.16153846
CMS2 0 214 2 4 0.02727273
CMS3 9 20 61 1 0.32967033
CMS4 4 9 0 162 0.07428571
# Vemos el summary del modelo summary(model_results$model) Vemos los resultados
# de la Matriz de Confusión de las predicciones (test)
print(model_results$results)
Confusion Matrix and Statistics
test_labels
predictions CMS1 CMS2 CMS3 CMS4
CMS1 58 0 1 2
CMS2 1 106 9 6
CMS3 2 1 35 0
CMS4 3 2 0 79
Overall Statistics
Accuracy : 0.9115
95% CI : (0.8738, 0.9409)
No Information Rate : 0.3574
P-Value [Acc > NIR] : < 2.2e-16
Kappa : 0.8767
Mcnemar's Test P-Value : NA
Statistics by Class:
Class: CMS1 Class: CMS2 Class: CMS3 Class: CMS4
Sensitivity 0.9062 0.9725 0.7778 0.9080
Specificity 0.9876 0.9184 0.9885 0.9771
Pos Pred Value 0.9508 0.8689 0.9211 0.9405
Neg Pred Value 0.9754 0.9836 0.9625 0.9638
Prevalence 0.2098 0.3574 0.1475 0.2852
Detection Rate 0.1902 0.3475 0.1148 0.2590
Detection Prevalence 0.2000 0.4000 0.1246 0.2754
Balanced Accuracy 0.9469 0.9454 0.8831 0.9426
# Regogemos los datos para comparar de la predicción
cms_Accuracy <- model_results$results$overall["Accuracy"]
cms_Kappa <- model_results$results$overall["Kappa"]
cms_AccuracyLower <- model_results$results$overall["AccuracyLower"]
cms_AccuracyUpper <- model_results$results$overall["AccuracyUpper"]
# Añadimos la variable 'class' que queramos para los modelos ML
df_exmat_cms_stage <- cbind(df_dataClin["cms_stage"], df_exmat_t)
# Cambiamos el nombre de la primera columna a 'class'
names(df_exmat_cms_stage)[1] <- "class"
# Ejecutamos el modelo Random Forest
model_results <- run_randomForest_model(df_exmat_cms_stage, ntree = 200)
# Vemos la información del modelo
model_results$model
Call:
randomForest(x = data_train[-1], y = train_labels, ntree = ntree)
Type of random forest: classification
Number of trees: 200
No. of variables tried at each split: 123
OOB estimate of error rate: 44.16%
Confusion matrix:
II_CMS1 II_CMS2 II_CMS3 II_CMS4 III_CMS1 III_CMS2 III_CMS3 III_CMS4
II_CMS1 66 1 4 8 3 1 0 0
II_CMS2 0 124 1 4 0 11 0 0
II_CMS3 9 11 43 0 0 1 0 0
II_CMS4 0 7 0 78 0 0 0 16
III_CMS1 28 5 1 5 7 0 0 1
III_CMS2 0 71 2 1 0 6 0 0
III_CMS3 1 6 16 2 0 0 2 0
III_CMS4 0 7 0 47 0 2 0 18
class.error
II_CMS1 0.2048193
II_CMS2 0.1142857
II_CMS3 0.3281250
II_CMS4 0.2277228
III_CMS1 0.8510638
III_CMS2 0.9250000
III_CMS3 0.9259259
III_CMS4 0.7567568
# Vemos el summary del modelo summary(model_results$model) Vemos los resultados
# de la Matriz de Confusión de las predicciones (test)
print(model_results$results)
Confusion Matrix and Statistics
test_labels
predictions II_CMS1 II_CMS2 II_CMS3 II_CMS4 III_CMS1 III_CMS2 III_CMS3 III_CMS4
II_CMS1 38 0 2 1 17 0 0 0
II_CMS2 0 66 8 2 2 35 2 1
II_CMS3 0 0 22 0 1 0 9 0
II_CMS4 2 2 0 42 2 0 0 30
III_CMS1 1 0 0 0 1 0 0 0
III_CMS2 0 2 0 0 0 4 2 0
III_CMS3 0 0 0 0 0 0 0 0
III_CMS4 0 0 0 5 0 0 0 6
Overall Statistics
Accuracy : 0.5869
95% CI : (0.5294, 0.6427)
No Information Rate : 0.2295
P-Value [Acc > NIR] : < 2.2e-16
Kappa : 0.4999
Mcnemar's Test P-Value : NA
Statistics by Class:
Class: II_CMS1 Class: II_CMS2 Class: II_CMS3
Sensitivity 0.9268 0.9429 0.68750
Specificity 0.9242 0.7872 0.96337
Pos Pred Value 0.6552 0.5690 0.68750
Neg Pred Value 0.9879 0.9788 0.96337
Prevalence 0.1344 0.2295 0.10492
Detection Rate 0.1246 0.2164 0.07213
Detection Prevalence 0.1902 0.3803 0.10492
Balanced Accuracy 0.9255 0.8650 0.82543
Class: II_CMS4 Class: III_CMS1 Class: III_CMS2
Sensitivity 0.8400 0.043478 0.10256
Specificity 0.8588 0.996454 0.98496
Pos Pred Value 0.5385 0.500000 0.50000
Neg Pred Value 0.9648 0.927393 0.88215
Prevalence 0.1639 0.075410 0.12787
Detection Rate 0.1377 0.003279 0.01311
Detection Prevalence 0.2557 0.006557 0.02623
Balanced Accuracy 0.8494 0.519966 0.54376
Class: III_CMS3 Class: III_CMS4
Sensitivity 0.00000 0.16216
Specificity 1.00000 0.98134
Pos Pred Value NaN 0.54545
Neg Pred Value 0.95738 0.89456
Prevalence 0.04262 0.12131
Detection Rate 0.00000 0.01967
Detection Prevalence 0.00000 0.03607
Balanced Accuracy 0.50000 0.57175
# Regogemos los datos para comparar de la predicción
cms_stage_Accuracy <- model_results$results$overall["Accuracy"]
cms_stage_Kappa <- model_results$results$overall["Kappa"]
cms_stage_AccuracyLower <- model_results$results$overall["AccuracyLower"]
cms_stage_AccuracyUpper <- model_results$results$overall["AccuracyUpper"]
# Filtramos y escogemos solo los datos msi_imputed = MSI
df_dataClin_exmat_msi <- df_dataClin_exmat_t[df_dataClin_exmat_t$msi_imputed == "MSI",
]
# Nos quedamos solo con la columna 'stage' como clase por lo que eliminamos las
# columnas de los datos clínicos que no nos interesan.
df_exmat_msi_stage <- df_dataClin_exmat_msi %>%
select(-ID, -msi_imputed, -cms, -grupo, -cms_stage, -study)
# Cambiamos el nombre de la primera columna a 'class'
names(df_exmat_msi_stage)[1] <- "class"
# Ejecutamos el modelo Random Forest
model_results <- run_randomForest_model(df_exmat_msi_stage)
# Vemos la información del modelo
model_results$model
Call:
randomForest(x = data_train[-1], y = train_labels, ntree = ntree)
Type of random forest: classification
Number of trees: 100
No. of variables tried at each split: 123
OOB estimate of error rate: 32.5%
Confusion matrix:
II III class.error
II 106 10 0.0862069
III 42 2 0.9545455
# Vemos el summary del modelo summary(model_results$model) Vemos los resultados
# de la Matriz de Confusión de las predicciones (test)
print(model_results$results)
Confusion Matrix and Statistics
test_labels
predictions II III
II 57 22
III 0 0
Accuracy : 0.7215
95% CI : (0.6093, 0.8165)
No Information Rate : 0.7215
P-Value [Acc > NIR] : 0.5571
Kappa : 0
Mcnemar's Test P-Value : 7.562e-06
Sensitivity : 1.0000
Specificity : 0.0000
Pos Pred Value : 0.7215
Neg Pred Value : NaN
Prevalence : 0.7215
Detection Rate : 0.7215
Detection Prevalence : 1.0000
Balanced Accuracy : 0.5000
'Positive' Class : II
# Regogemos los datos para comparar de la predicción
msi_stage_Accuracy <- model_results$results$overall["Accuracy"]
msi_stage_Kappa <- model_results$results$overall["Kappa"]
msi_stage_AccuracyLower <- model_results$results$overall["AccuracyLower"]
msi_stage_AccuracyUpper <- model_results$results$overall["AccuracyUpper"]
# Filtramos y escogemos solo los datos msi_imputed = MSS
df_dataClin_exmat_mss <- df_dataClin_exmat_t[df_dataClin_exmat_t$msi_imputed == "MSS",
]
# Nos quedamos solo con la columna 'stage' como clase por lo que eliminamos las
# columnas de los datos clínicos que no nos interesan.
df_exmat_mss_stage <- df_dataClin_exmat_mss %>%
select(-ID, -msi_imputed, -cms, -grupo, -cms_stage, -study)
# Cambiamos el nombre de la primera columna a 'class'
names(df_exmat_mss_stage)[1] <- "class"
# Ejecutamos el modelo Random Forest
model_results <- run_randomForest_model(df_exmat_mss_stage)
# Vemos la información del modelo
model_results$model
Call:
randomForest(x = data_train[-1], y = train_labels, ntree = ntree)
Type of random forest: classification
Number of trees: 100
No. of variables tried at each split: 123
OOB estimate of error rate: 39.57%
Confusion matrix:
II III class.error
II 287 54 0.1583578
III 168 52 0.7636364
# Vemos el summary del modelo summary(model_results$model) Vemos los resultados
# de la Matriz de Confusión de las predicciones (test)
print(model_results$results)
Confusion Matrix and Statistics
test_labels
predictions II III
II 155 86
III 15 23
Accuracy : 0.638
95% CI : (0.5786, 0.6944)
No Information Rate : 0.6093
P-Value [Acc > NIR] : 0.1789
Kappa : 0.139
Mcnemar's Test P-Value : 3.278e-12
Sensitivity : 0.9118
Specificity : 0.2110
Pos Pred Value : 0.6432
Neg Pred Value : 0.6053
Prevalence : 0.6093
Detection Rate : 0.5556
Detection Prevalence : 0.8638
Balanced Accuracy : 0.5614
'Positive' Class : II
# Regogemos los datos para comparar de la predicción
mss_stage_Accuracy <- model_results$results$overall["Accuracy"]
mss_stage_Kappa <- model_results$results$overall["Kappa"]
mss_stage_AccuracyLower <- model_results$results$overall["AccuracyLower"]
mss_stage_AccuracyUpper <- model_results$results$overall["AccuracyUpper"]
# Filtramos y escogemos solo los datos msi_imputed = MSS
df_dataClin_exmat_cms1 <- df_dataClin_exmat_t[df_dataClin_exmat_t$cms == "CMS1",
]
# Nos quedamos solo con la columna 'stage' como clase por lo que eliminamos las
# columnas de los datos clínicos que no nos interesan.
df_exmat_cms1_stage <- df_dataClin_exmat_cms1 %>%
select(-ID, -msi_imputed, -cms, -grupo, -cms_stage, -study)
# Cambiamos el nombre de la primera columna a 'class'
names(df_exmat_cms1_stage)[1] <- "class"
# Ejecutamos el modelo Random Forest
model_results <- run_randomForest_model(df_exmat_cms1_stage)
# Vemos la información del modelo
model_results$model
Call:
randomForest(x = data_train[-1], y = train_labels, ntree = ntree)
Type of random forest: classification
Number of trees: 100
No. of variables tried at each split: 123
OOB estimate of error rate: 36.15%
Confusion matrix:
II III class.error
II 73 10 0.1204819
III 37 10 0.7872340
# Vemos el summary del modelo summary(model_results$model) Vemos los resultados
# de la Matriz de Confusión de las predicciones (test)
print(model_results$results)
Confusion Matrix and Statistics
test_labels
predictions II III
II 39 19
III 2 4
Accuracy : 0.6719
95% CI : (0.5431, 0.7841)
No Information Rate : 0.6406
P-Value [Acc > NIR] : 0.3520305
Kappa : 0.1494
Mcnemar's Test P-Value : 0.0004803
Sensitivity : 0.9512
Specificity : 0.1739
Pos Pred Value : 0.6724
Neg Pred Value : 0.6667
Prevalence : 0.6406
Detection Rate : 0.6094
Detection Prevalence : 0.9062
Balanced Accuracy : 0.5626
'Positive' Class : II
# Regogemos los datos para comparar de la predicción
cms1_stage_Accuracy <- model_results$results$overall["Accuracy"]
cms1_stage_Kappa <- model_results$results$overall["Kappa"]
cms1_stage_AccuracyLower <- model_results$results$overall["AccuracyLower"]
cms1_stage_AccuracyUpper <- model_results$results$overall["AccuracyUpper"]
# Filtramos y escogemos solo los datos msi_imputed = MSS
df_dataClin_exmat_cms2 <- df_dataClin_exmat_t[df_dataClin_exmat_t$cms == "CMS2",
]
# Nos quedamos solo con la columna 'stage' como clase por lo que eliminamos las
# columnas de los datos clínicos que no nos interesan.
df_exmat_cms2_stage <- df_dataClin_exmat_cms2 %>%
select(-ID, -msi_imputed, -cms, -grupo, -cms_stage, -study)
# Cambiamos el nombre de la primera columna a 'class'
names(df_exmat_cms2_stage)[1] <- "class"
# Ejecutamos el modelo Random Forest
model_results <- run_randomForest_model(df_exmat_cms2_stage)
# Vemos la información del modelo
model_results$model
Call:
randomForest(x = data_train[-1], y = train_labels, ntree = ntree)
Type of random forest: classification
Number of trees: 100
No. of variables tried at each split: 123
OOB estimate of error rate: 40.45%
Confusion matrix:
II III class.error
II 124 16 0.1142857
III 73 7 0.9125000
# Vemos el summary del modelo summary(model_results$model) Vemos los resultados
# de la Matriz de Confusión de las predicciones (test)
print(model_results$results)
Confusion Matrix and Statistics
test_labels
predictions II III
II 65 36
III 5 3
Accuracy : 0.6239
95% CI : (0.526, 0.7148)
No Information Rate : 0.6422
P-Value [Acc > NIR] : 0.6937
Kappa : 0.0067
Mcnemar's Test P-Value : 2.797e-06
Sensitivity : 0.92857
Specificity : 0.07692
Pos Pred Value : 0.64356
Neg Pred Value : 0.37500
Prevalence : 0.64220
Detection Rate : 0.59633
Detection Prevalence : 0.92661
Balanced Accuracy : 0.50275
'Positive' Class : II
# Regogemos los datos para comparar de la predicción
cms2_stage_Accuracy <- model_results$results$overall["Accuracy"]
cms2_stage_Kappa <- model_results$results$overall["Kappa"]
cms2_stage_AccuracyLower <- model_results$results$overall["AccuracyLower"]
cms2_stage_AccuracyUpper <- model_results$results$overall["AccuracyUpper"]
# Filtramos y escogemos solo los datos msi_imputed = MSS
df_dataClin_exmat_cms3 <- df_dataClin_exmat_t[df_dataClin_exmat_t$cms == "CMS3",
]
# Nos quedamos solo con la columna 'stage' como clase por lo que eliminamos las
# columnas de los datos clínicos que no nos interesan.
df_exmat_cms3_stage <- df_dataClin_exmat_cms3 %>%
select(-ID, -msi_imputed, -cms, -grupo, -cms_stage, -study)
# Cambiamos el nombre de la primera columna a 'class'
names(df_exmat_cms3_stage)[1] <- "class"
# Ejecutamos el modelo Random Forest
model_results <- run_randomForest_model(df_exmat_cms3_stage)
# Vemos la información del modelo
model_results$model
Call:
randomForest(x = data_train[-1], y = train_labels, ntree = ntree)
Type of random forest: classification
Number of trees: 100
No. of variables tried at each split: 123
OOB estimate of error rate: 28.57%
Confusion matrix:
II III class.error
II 63 1 0.0156250
III 25 2 0.9259259
# Vemos el summary del modelo summary(model_results$model) Vemos los resultados
# de la Matriz de Confusión de las predicciones (test)
print(model_results$results)
Confusion Matrix and Statistics
test_labels
predictions II III
II 31 13
III 1 0
Accuracy : 0.6889
95% CI : (0.5335, 0.8183)
No Information Rate : 0.7111
P-Value [Acc > NIR] : 0.695081
Kappa : -0.043
Mcnemar's Test P-Value : 0.003283
Sensitivity : 0.9688
Specificity : 0.0000
Pos Pred Value : 0.7045
Neg Pred Value : 0.0000
Prevalence : 0.7111
Detection Rate : 0.6889
Detection Prevalence : 0.9778
Balanced Accuracy : 0.4844
'Positive' Class : II
# Regogemos los datos para comparar de la predicción
cms3_stage_Accuracy <- model_results$results$overall["Accuracy"]
cms3_stage_Kappa <- model_results$results$overall["Kappa"]
cms3_stage_AccuracyLower <- model_results$results$overall["AccuracyLower"]
cms3_stage_AccuracyUpper <- model_results$results$overall["AccuracyUpper"]
# Filtramos y escogemos solo los datos msi_imputed = MSS
df_dataClin_exmat_cms4 <- df_dataClin_exmat_t[df_dataClin_exmat_t$cms == "CMS4",
]
# Nos quedamos solo con la columna 'stage' como clase por lo que eliminamos las
# columnas de los datos clínicos que no nos interesan.
df_exmat_cms4_stage <- df_dataClin_exmat_cms4 %>%
select(-ID, -msi_imputed, -cms, -grupo, -cms_stage, -study)
# Cambiamos el nombre de la primera columna a 'class'
names(df_exmat_cms4_stage)[1] <- "class"
# Ejecutamos el modelo Random Forest
model_results <- run_randomForest_model(df_exmat_cms4_stage)
# Vemos la información del modelo
model_results$model
Call:
randomForest(x = data_train[-1], y = train_labels, ntree = ntree)
Type of random forest: classification
Number of trees: 100
No. of variables tried at each split: 123
OOB estimate of error rate: 42.86%
Confusion matrix:
II III class.error
II 76 25 0.2475248
III 50 24 0.6756757
# Vemos el summary del modelo summary(model_results$model) Vemos los resultados
# de la Matriz de Confusión de las predicciones (test)
print(model_results$results)
Confusion Matrix and Statistics
test_labels
predictions II III
II 42 30
III 8 7
Accuracy : 0.5632
95% CI : (0.4526, 0.6694)
No Information Rate : 0.5747
P-Value [Acc > NIR] : 0.6292447
Kappa : 0.0316
Mcnemar's Test P-Value : 0.0006577
Sensitivity : 0.8400
Specificity : 0.1892
Pos Pred Value : 0.5833
Neg Pred Value : 0.4667
Prevalence : 0.5747
Detection Rate : 0.4828
Detection Prevalence : 0.8276
Balanced Accuracy : 0.5146
'Positive' Class : II
# Regogemos los datos para comparar de la predicción
cms4_stage_Accuracy <- model_results$results$overall["Accuracy"]
cms4_stage_Kappa <- model_results$results$overall["Kappa"]
cms4_stage_AccuracyLower <- model_results$results$overall["AccuracyLower"]
cms4_stage_AccuracyUpper <- model_results$results$overall["AccuracyUpper"]
| Modelo | Accuraccy | kappa | 95% CI |
|---|---|---|---|
| Estadio | 0.6462396 | 0.0870262 | (0.5943389 - 0.6957054) |
| MSI/MSS | 0.8551532 | 0.4899732 | (0.8144282 - 0.8898921) |
| MSI/MSS-Estadio | 0.5530726 | 0.2427587 | (0.4999164 - 0.6053439) |
| CMS | 0.9114754 | 0.8767308 | (0.8738159 - 0.9408505) |
| CMS-Estadio | 0.5868852 | 0.4998568 | (0.5293643 - 0.6427016) |
| MSI-Estadio | 0.721519 | 0 | (0.6092649 - 0.816545) |
| MSS-Estadio | 0.6379928 | 0.1390204 | (0.5785781 - 0.6944416) |
| CMS1-Estadio | 0.671875 | 0.1493671 | (0.5431232 - 0.784128) |
| CMS2-Estadio | 0.6238532 | 0.0066681 | (0.525959 - 0.7148372) |
| CMS3-Estadio | 0.6888889 | -0.0430464 | (0.533509 - 0.8183412) |
| CMS4-Estadio | 0.5632184 | 0.0316344 | (0.4526478 - 0.669361) |
Repositorio: (Ferrón 2024)